# 9/19  run decomp models using specification from 230918


library(RSiena)
library(parallel)

#Set working directory
setwd('/Users/drschaef/Google Drive/Projects/Adriana/Papers/Homophily/Output')
setwd('/Users/drschaef/My Drive/Projects/Adriana/Papers/Homophily/Output')

rm(list=ls())

mydatRF <-  readRDS("rfnetM123_221014.rds") # friendships and race only

# effect-setting functions
set0_raceMain <- function(EFF_OBJ) {
	# main effects of race - excluded category in White
	EFF_OBJ <- includeEffects(EFF_OBJ, altX, name="raceNet", interaction1="black") 
	EFF_OBJ <- includeEffects(EFF_OBJ, altX, name="raceNet", interaction1="natam") 
	EFF_OBJ <- includeEffects(EFF_OBJ, altX, name="raceNet", interaction1="latin")
	EFF_OBJ <- includeEffects(EFF_OBJ, altX, name="raceNet", interaction1="asian") 
	EFF_OBJ <- includeEffects(EFF_OBJ, altX, name="raceNet", interaction1="orace")
	EFF_OBJ <- includeEffects(EFF_OBJ, altX, name="raceNet", interaction1="mrace") 
	}
set1_raceControls_alter <- function(EFF_OBJ, ATTR, INCL) {
	EFF_OBJ <- includeEffects(EFF_OBJ, egoX, name="raceNet", interaction1=ATTR)
	EFF_OBJ <- includeInteraction(EFF_OBJ, egoX, altX, name="raceNet", interaction1=c(ATTR,'mrace'), include=INCL)
	EFF_OBJ <- includeInteraction(EFF_OBJ, egoX, altX, name="raceNet", interaction1=c(ATTR,'asian'), include=INCL)
	EFF_OBJ <- includeInteraction(EFF_OBJ, egoX, altX, name="raceNet", interaction1=c(ATTR,'black'), include=INCL)
	EFF_OBJ <- includeInteraction(EFF_OBJ, egoX, altX, name="raceNet", interaction1=c(ATTR,'latin'), include=INCL)
	EFF_OBJ <- includeInteraction(EFF_OBJ, egoX, altX, name="raceNet", interaction1=c(ATTR,'orace'), include=INCL)
	}
set1_raceControls_simDist2 <- function(EFF_OBJ, ATTR, INCL) {
	EFF_OBJ <- includeEffects(EFF_OBJ, egoX, name="raceNet", interaction1=ATTR, include=INCL)
	EFF_OBJ <- includeEffects(EFF_OBJ, simEgoInDist2, name="raceNet", interaction1=ATTR, include=INCL)
	}
set1_raceControls_selection <- function(EFF_OBJ, INCL) {
	EFF_OBJ <- setEffect(EFF_OBJ, sameX, name="friendship", interaction1="female", include=INCL, initialValue=.47)
	EFF_OBJ <- setEffect(EFF_OBJ, sameX, name="friendship", interaction1="gradecohort", include=INCL, initialValue=1.06)
	EFF_OBJ <- setEffect(EFF_OBJ, simX, name="friendship", interaction1="pared8", include=INCL, initialValue=.21)
	EFF_OBJ <- setEffect(EFF_OBJ, sameX, name="friendship", interaction1="immgen", include=INCL, initialValue=.13)
	EFF_OBJ <- setEffect(EFF_OBJ, simX, name="friendship", interaction1="acad9", include=INCL, initialValue=.5)
#	EFF_OBJ <- setEffect(EFF_OBJ, simX, name="friendship", interaction1="perla", include=INCL, initialValue=.3)
	EFF_OBJ <- setEffect(EFF_OBJ, simX, name="friendship", interaction1="disc", include=INCL, initialValue=.16)
	}
set2_ecaEffects <- function(EFF_OBJ, INCL) {
	# restrict eca co-evolution network effects to 0 and fixed
	EFF_OBJ$fix[EFF_OBJ$include==T & EFF_OBJ$name=='ecaEndog'] <- TRUE
	EFF_OBJ$initialValue[EFF_OBJ$include==T & EFF_OBJ$name=='ecaEndog'] <- 0
	EFF_OBJ$initialValue[EFF_OBJ$include==T & EFF_OBJ$name=='ecaEndog'][1:2] <- .0001  # rate effects can't be 0
	# eca influence
	EFF_OBJ <- setEffect(EFF_OBJ, to, name="raceNet", interaction1="ecaEndog", incl=INCL)
	}
set2_ecaFriendEffects <- function(EFF_OBJ, INCL) {
	# eca selection
	EFF_OBJ <- setEffect(EFF_OBJ, X, name="friendship", interaction1="ecaDum", include=INCL, initialValue=.63)
	}	
set3_raceInfluence <- function(EFF_OBJ) {
	EFF_OBJ <- setEffect(EFF_OBJ, to, name="raceNet", interaction1="friendship", initialValue=.6)  # friend influence
	EFF_OBJ <- setEffect(EFF_OBJ, to, name="raceNet", interaction1="friendship", type='creation', initialValue=-.3)  # friend influence
	}
set4_ancestry <- function(EFF_OBJ, INCL) {
	EFF_OBJ <- setEffect(EFF_OBJ, X, name="friendship", interaction1="mumu", include=INCL, initialValue=.3)
	EFF_OBJ <- setEffect(EFF_OBJ, X, name="friendship", interaction1="mumo", include=INCL, initialValue=.08)
	EFF_OBJ <- setEffect(EFF_OBJ, X, name="friendship", interaction1="momu", include=INCL, initialValue=.16)
	}
set5_endogNet <- function(EFF_OBJ, INCL) {
	EFF_OBJ <- setEffect(EFF_OBJ, recip, name="friendship", include=INCL, initialValue=3.3)
	EFF_OBJ <- setEffect(EFF_OBJ, outActSqrt, name="friendship", include=INCL, initialValue=.29)
	EFF_OBJ <- setEffect(EFF_OBJ, inPopSqrt, name="friendship", include=INCL, initialValue=.4) 
	EFF_OBJ <- setEffect(EFF_OBJ, outPopSqrt, name="friendship", include=INCL, initialValue=-.63)
	EFF_OBJ <- setEffect(EFF_OBJ, transTrip, name='friendship', include=INCL, initialValue=.76)
	EFF_OBJ <- setEffect(EFF_OBJ, nbrDist2, name="friendship", include=INCL, initialValue=-.09)
	EFF_OBJ <- setEffect(EFF_OBJ, transRecTrip, name="friendship", include=INCL, initialValue=-.76)
	}
set6_raceSelection <- function(EFF_OBJ) {
	# choose friends with same self-classification
	EFF_OBJ <- setEffect(EFF_OBJ, from, name="friendship", interaction1="raceNet", include=T, initialValue=.7)
	EFF_OBJ <- setEffect(EFF_OBJ, from, name="friendship", interaction1="raceNet", type='endow', include=T, initialValue=-.1)
	}

myeff_trend <- getEffects( mydatRF, nintn=60 )   
myeff_trend <- setEffect(myeff_trend, recip, name='friendship', include=F)

# specify common model effects
myeff0 <- getEffects( mydatRF, nintn=60 )   # race and friendship network
myeff0 <- set0_raceMain(myeff0)
myeff0 <- setEffect(myeff0, RateX, type='rate', name="raceNet", interaction1="mono", fix=T, initialValue=-100) 
myeff0 <- includeTimeDummy(myeff0, density, timeDummy = "all")
myeff0$initialValue[myeff0$include==T][1:4] <- c(9, 8.5, -5.5, 3.3)

myeff_control <- myeff0
myeff_control <- set1_raceControls_alter(myeff_control, 'perla', T)
myeff_control <- set1_raceControls_alter(myeff_control, 'disc', T)
myeff_control <- set1_raceControls_simDist2(myeff_control, 'perla', T)
myeff_control <- set1_raceControls_simDist2(myeff_control, 'disc', T)
myeff_control <- set1_raceControls_selection(myeff_control, T)

myeff_control <- set2_ecaFriendEffects(myeff_control, T)
#myeff_control <- set4_ancestry(myeff_control, T)
myeff_control <- set5_endogNet(myeff_control, T)

myeff_full <- myeff_control
#myeff_full <- set3_raceInfluence(myeff_full)
#myeff_full <- set6_raceSelection(myeff_full)
myeff_full <- setEffect(myeff_full, from, name="friendship", interaction1="raceNet", type='eval', include=F, fix=F, test=F)
myeff_full <- setEffect(myeff_full, from, name="friendship", interaction1="raceNet", type='creation', include=T, fix=F, test=F)
myeff_full <- setEffect(myeff_full, from, name="friendship", interaction1="raceNet", type='endow', include=T, fix=F, test=F)
myeff_full <- setEffect(myeff_full, X, name="friendship", interaction1="mumu", type='eval', include=F, fix=F, test=F)
myeff_full <- setEffect(myeff_full, X, name="friendship", interaction1="mumo", type='eval', include=F, fix=F, test=F)
myeff_full <- setEffect(myeff_full, X, name="friendship", interaction1="momu", type='eval', include=F, fix=F, test=F)
myeff_full <- setEffect(myeff_full, X, name="friendship", interaction1="mumu", type='creation', include=T, fix=F, test=F)
myeff_full <- setEffect(myeff_full, X, name="friendship", interaction1="mumo", type='creation', include=T, fix=F, test=F)
myeff_full <- setEffect(myeff_full, X, name="friendship", interaction1="momu", type='creation', include=T, fix=F, test=F)
myeff_full <- setEffect(myeff_full, X, name="friendship", interaction1="mumu", type='endow', include=T, fix=F, test=F)
myeff_full <- setEffect(myeff_full, X, name="friendship", interaction1="mumo", type='endow', include=T, fix=F, test=F)
myeff_full <- setEffect(myeff_full, X, name="friendship", interaction1="momu", type='endow', include=T, fix=F, test=F)
myeff_full <- setEffect(myeff_full, to, name="raceNet", interaction1="friendship", type='eval', include=F, fix=F, test=F)
myeff_full <- setEffect(myeff_full, to, name="raceNet", interaction1="friendship", type='creation', include=T, fix=F, test=F)
myeff_full <- setEffect(myeff_full, to, name="raceNet", interaction1="friendship", type='endow', include=T, fix=F, test=F)

##define the algorithm settings
setwd('/Users/drschaef/Google Drive/Projects/Adriana/Papers/Homophily/Output/tempM')
setwd('/Users/drschaef/My Drive/Projects/Adriana/Papers/Homophily/Output/tempM')
myAlgorithm <- sienaAlgorithmCreate( projname = 'model230927_decomp_M', nsub=4, n3=2000)

# get starting values
m_full <- readRDS('fullModel_M_230927.RDS')  # load in the full model 

# fitted models needed
m_trend <- siena07( myAlgorithm, data = mydatRF, effects = myeff_trend, useCluster=T, initC=T, nbrNodes=7, prevAns=m_full, returnDeps=T)
m_control <- siena07( myAlgorithm, data = mydatRF, effects = myeff_control, useCluster=T, initC=T, nbrNodes=7, prevAns=m_full, returnDeps=T)
saveRDS(m_trend, 'model230927_trend_M.RDS')
saveRDS(m_control, 'model230927_control_M.RDS')

m_trend <- readRDS('model230927_trend_M.RDS')  # load in the full model (which is technically reduced)
m_control <- readRDS('model230927_control_M.RDS')  # load in the full model (which is technically reduced)

# simulate fitted and hybrid
# - fitted full (reduced) model
simAlgorithm <- sienaAlgorithmCreate( projname = 'model230927_sim_M', cond=F, useStdInits=F, nsub=0, simOnly=T, n3=1000)
myeff_fullUpdated <- updateTheta(myeff_full, m_full)
s_full <- siena07( simAlgorithm, data = mydatRF, effects = myeff_fullUpdated, useCluster=T, initC=T, nbrNodes=7, returnDeps=T)
saveRDS(s_full, 'model230927_simFull_M.RDS')
s_full <- readRDS('model230927_simFull_M.RDS')  # load in the full model (which is technically reduced)

# - hybrid models
#  - influence all, selection creation
myeff_inflAll_selCre <- myeff_full   
myeff_inflAll_selCre <- setEffect(myeff_inflAll_selCre, from, name="friendship", interaction1="raceNet", type='endow', include=F, fix=F, test=F)
myeff_inflAll_selCre <- setEffect(myeff_inflAll_selCre, X, name="friendship", interaction1="mumu", type='endow', include=F, fix=F, test=F)
myeff_inflAll_selCre <- setEffect(myeff_inflAll_selCre, X, name="friendship", interaction1="mumo", type='endow', include=F, fix=F, test=F)
myeff_inflAll_selCre <- setEffect(myeff_inflAll_selCre, X, name="friendship", interaction1="momu", type='endow', include=F, fix=F, test=F)
myeff_inflAll_selCre <- updateTheta(myeff_inflAll_selCre, m_full)   

#  - influence all, selection endowment
myeff_inflAll_selEnd <- myeff_full   
myeff_inflAll_selEnd <- setEffect(myeff_inflAll_selEnd, from, name="friendship", interaction1="raceNet", type='creation', include=F, fix=F, test=F)
myeff_inflAll_selEnd <- setEffect(myeff_inflAll_selEnd, X, name="friendship", interaction1="mumu", type='creation', include=F, fix=F, test=F)
myeff_inflAll_selEnd <- setEffect(myeff_inflAll_selEnd, X, name="friendship", interaction1="mumo", type='creation', include=F, fix=F, test=F)
myeff_inflAll_selEnd <- setEffect(myeff_inflAll_selEnd, X, name="friendship", interaction1="momu", type='creation', include=F, fix=F, test=F)
myeff_inflAll_selEnd <- updateTheta(myeff_inflAll_selEnd, m_full)   

#  - influence all, selection none
myeff_inflAll_selNon <- myeff_full   
myeff_inflAll_selNon <- setEffect(myeff_inflAll_selNon, from, name="friendship", interaction1="raceNet", type='creation', include=F, fix=F, test=F)
myeff_inflAll_selNon <- setEffect(myeff_inflAll_selNon, X, name="friendship", interaction1="mumu", type='creation', include=F, fix=F, test=F)
myeff_inflAll_selNon <- setEffect(myeff_inflAll_selNon, X, name="friendship", interaction1="mumo", type='creation', include=F, fix=F, test=F)
myeff_inflAll_selNon <- setEffect(myeff_inflAll_selNon, X, name="friendship", interaction1="momu", type='creation', include=F, fix=F, test=F)
myeff_inflAll_selNon <- setEffect(myeff_inflAll_selNon, from, name="friendship", interaction1="raceNet", type='endow', include=F, fix=F, test=F)
myeff_inflAll_selNon <- setEffect(myeff_inflAll_selNon, X, name="friendship", interaction1="mumu", type='endow', include=F, fix=F, test=F)
myeff_inflAll_selNon <- setEffect(myeff_inflAll_selNon, X, name="friendship", interaction1="mumo", type='endow', include=F, fix=F, test=F)
myeff_inflAll_selNon <- setEffect(myeff_inflAll_selNon, X, name="friendship", interaction1="momu", type='endow', include=F, fix=F, test=F)
myeff_inflAll_selNon <- updateTheta(myeff_inflAll_selNon, m_full)   

#  - influence creation, selection all
myeff_inflCre_selAll <- myeff_full   
myeff_inflCre_selAll <- setEffect(myeff_inflCre_selAll, to, name="raceNet", interaction1="friendship", type='endow', include=F, fix=F, test=F)
myeff_inflCre_selAll <- updateTheta(myeff_inflCre_selAll, m_full)   

#  - influence endowment, selection all
myeff_inflEnd_selAll <- myeff_full   
myeff_inflEnd_selAll <- setEffect(myeff_inflEnd_selAll, to, name="raceNet", interaction1="friendship", type='creation', include=F, fix=F, test=F)
myeff_inflEnd_selAll <- updateTheta(myeff_inflEnd_selAll, m_full)   

#  - influence none, selection all
myeff_inflNon_selAll <- myeff_full   
myeff_inflNon_selAll <- setEffect(myeff_inflNon_selAll, to, name="raceNet", interaction1="friendship", type='creation', include=F, fix=F, test=F)
myeff_inflNon_selAll <- setEffect(myeff_inflNon_selAll, to, name="raceNet", interaction1="friendship", type='endow', include=F, fix=F, test=F)
myeff_inflNon_selAll <- updateTheta(myeff_inflNon_selAll, m_full)   

#  - simulate all 6 models
s_inflAll_selCre <- siena07( simAlgorithm, data = mydatRF, effects = myeff_inflAll_selCre, useCluster=T, initC=T, nbrNodes=7, returnDeps=T)
s_inflAll_selEnd <- siena07( simAlgorithm, data = mydatRF, effects = myeff_inflAll_selEnd, useCluster=T, initC=T, nbrNodes=7, returnDeps=T)
s_inflCre_selAll <- siena07( simAlgorithm, data = mydatRF, effects = myeff_inflCre_selAll, useCluster=T, initC=T, nbrNodes=7, returnDeps=T)
s_inflEnd_selAll <- siena07( simAlgorithm, data = mydatRF, effects = myeff_inflEnd_selAll, useCluster=T, initC=T, nbrNodes=7, returnDeps=T)
s_inflAll_selNon <- siena07( simAlgorithm, data = mydatRF, effects = myeff_inflAll_selNon, useCluster=T, initC=T, nbrNodes=7, returnDeps=T)
s_inflNon_selAll <- siena07( simAlgorithm, data = mydatRF, effects = myeff_inflNon_selAll, useCluster=T, initC=T, nbrNodes=7, returnDeps=T)

save.image('/Users/drschaef/Desktop/workspaceM230927.RData')
load('/Users/drschaef/Desktop/workspaceM230927.RData')

# sim data structure
# $sims[[a]][[b]][[c]][[d]]
# a = list, each element is one simulation run (2002) 
# b = list, each element is one period (e.g., Data1, Data2)
# c = list, each element is a dependent network (e.g., friendship, raceNet, ecaEndog)
# d = list, one element, named "1" and contains an edgelist
# begin the decomp
require(network)
require(sna)
require(netseg)
library(Matrix)

# a toy model for testing
m_trend <- readRDS('/Users/drschaef/Google Drive/Projects/Adriana/Papers/Homophily/Output/tempM/model230927_trend_M.RDS')
myeff_trendUpdated <- updateTheta(myeff_trend, m_trend)
simAlgorithm30 <- sienaAlgorithmCreate( projname = 'temp', cond=F, useStdInits=F, nsub=0, simOnly=T, n3=30)
s_trend <- siena07( simAlgorithm30, data = mydatRF, effects = myeff_trendUpdated, returnDeps=T)

# recode race into single measure; multiple race get dropped
# but first, count how many have 0 race, more than 1 race
#lapply(simrace20, function(x) table(degree(x, cmode='outdegree')))
extractFrNet <- function(s, DAT) {
	nS <- length(DAT$Data1$nodeSets$students)
	net <- network.initialize(n=nS, directed=TRUE)
	edges <- s$'Data1'$'friendship'$'1'[,1:2]
	net <- network.edgelist(edges,net)
  	return(net)
	}
extractRaNet <- function(s, DAT) {
	nS <- length(DAT$Data1$nodeSets$students)
	nR <- length(DAT$Data1$nodeSets$erGroups)
	edges <- s$'Data1'$'raceNet'$'1'  
	netDims <- c(nS, nR)  # two mode network dimensions
	edgesPlus <- rbind(edges,c(netDims,0))    # add a final row that ensures the correct network dimensions
	A <- sparseMatrix(i=edgesPlus[,1], j= edgesPlus[,2], x= edgesPlus[,3]) 
	mat <- as.matrix(A)
	return(mat)
	}
calcKSEI <- function(SIMNET, DAT) {
	simnets <- lapply(SIMNET$sims, function(x) extractFrNet(x,mydatRF))
	simrace <- lapply(SIMNET$sims, function(x) extractRaNet(x,mydatRF))
	KSEI <- mapply(x=simnets, y=simrace, function(x,y)  {
		nS <- length(DAT$Data1$nodeSets$students)
		od <- degree(y, cmode='outdegree')[1:nS]   # only keep cases reporting one race
		keep <- od==1
		netKeep <- x[keep==1,keep==1]
		raceKeep <- y[keep==1,]
		raceCat <- apply(raceKeep, 1, function(z) which(z==1))
		net <- network(netKeep)
		set.vertex.attribute(net,'race',raceCat)
		mm <- mixingmatrix(net,'race')
		ksei <- ei(mm)
		return(ksei)		
		})
	return(KSEI) 
	}

ksei_trend <- calcKSEI(s_trend, mydatRF)  # test
ksei_trend <- calcKSEI(m_trend, mydatRF)  # fitted model 
ksei_control <- calcKSEI(m_control, mydatRF)  # fitted model
ksei_full <- calcKSEI(s_full, mydatRF)  # first object is full answer object
ksei_inflAll_selCre <- calcKSEI(s_inflAll_selCre, mydatRF)  # first object is full answer object
ksei_inflAll_selEnd <- calcKSEI(s_inflAll_selEnd, mydatRF)  # first object is full answer object
ksei_inflAll_selNon <- calcKSEI(s_inflAll_selNon, mydatRF)  # first object is full answer object
ksei_inflCre_selAll <- calcKSEI(s_inflCre_selAll, mydatRF)  # first object is full answer object
ksei_inflEnd_selAll <- calcKSEI(s_inflEnd_selAll, mydatRF)  # first object is full answer object
ksei_inflNon_selAll <- calcKSEI(s_inflNon_selAll, mydatRF)  # first object is full answer object

setwd('/Users/drschaef/Google Drive/Projects/Adriana/Papers/Homophily/Output')
setwd('/Users/drschaef/My Drive/Projects/Adriana/Papers/Homophily/Output')
#save(list=ls()[ls() %in% ls(pattern="^ksei")], file='ksei_M230918.RData')

# REDO WITH SECOND PERIOD (change to Data2)
# recode race into single measure; multiple race get dropped
# but first, count how many have 0 race, more than 1 race
#lapply(simrace20, function(x) table(degree(x, cmode='outdegree')))
extractFrNet_d2 <- function(s, DAT) {
	nS <- length(DAT$Data2$nodeSets$students)
	net <- network.initialize(n=nS, directed=TRUE)
	edges <- s$'Data2'$'friendship'$'3'[,1:2]
	net <- network.edgelist(edges,net)
  	return(net)
	}
extractRaNet_d2 <- function(s, DAT) {
	nS <- length(DAT$Data2$nodeSets$students)
	nR <- length(DAT$Data2$nodeSets$erGroups)
	edges <- s$'Data2'$'raceNet'$'3'  
	netDims <- c(nS, nR)  # two mode network dimensions
	edgesPlus <- rbind(edges,c(netDims,0))    # add a final row that ensures the correct network dimensions
	A <- sparseMatrix(i=edgesPlus[,1], j= edgesPlus[,2], x= edgesPlus[,3]) 
	mat <- as.matrix(A)
	return(mat)
	}
calcKSEI_d2 <- function(SIMNET, DAT) {
	simnets <- lapply(SIMNET$sims, function(x) extractFrNet_d2(x,mydatRF))
	simrace <- lapply(SIMNET$sims, function(x) extractRaNet_d2(x,mydatRF))
	KSEI <- mapply(x=simnets, y=simrace, function(x,y)  {
		nS <- length(DAT$Data2$nodeSets$students)
		od <- degree(y, cmode='outdegree')[1:nS]   # only keep cases reporting one race
		keep <- od==1
		netKeep <- x[keep==1,keep==1]
		raceKeep <- y[keep==1,]
		raceCat <- apply(raceKeep, 1, function(z) which(z==1))
		net <- network(netKeep)
		set.vertex.attribute(net,'race',raceCat)
		mm <- mixingmatrix(net,'race')
		ksei <- ei(mm)
		return(ksei)		
		})
	return(KSEI) 
	}

#ksei_trend <- calcKSEI(s_trend, mydatRF)  # test
ksei_trend_d2 <- calcKSEI_d2(m_trend, mydatRF)  # fitted model 
ksei_control_d2 <- calcKSEI_d2(m_control, mydatRF)  # fitted model
ksei_full_d2 <- calcKSEI_d2(s_full, mydatRF)  # first object is full answer object
ksei_inflAll_selCre_d2 <- calcKSEI_d2(s_inflAll_selCre, mydatRF)  # first object is full answer object
ksei_inflAll_selEnd_d2 <- calcKSEI_d2(s_inflAll_selEnd, mydatRF)  # first object is full answer object
ksei_inflAll_selNon_d2 <- calcKSEI_d2(s_inflAll_selNon, mydatRF)  # first object is full answer object
ksei_inflCre_selAll_d2 <- calcKSEI_d2(s_inflCre_selAll, mydatRF)  # first object is full answer object
ksei_inflEnd_selAll_d2 <- calcKSEI_d2(s_inflEnd_selAll, mydatRF)  # first object is full answer object
ksei_inflNon_selAll_d2 <- calcKSEI_d2(s_inflNon_selAll, mydatRF)  # first object is full answer object

setwd('/Users/drschaef/Google Drive/Projects/Adriana/Papers/Homophily/Output')
setwd('/Users/drschaef/My Drive/Projects/Adriana/Papers/Homophily/Output')
#save(list=ls()[ls() %in% ls(pattern="^ksei")], file='ksei_M230918_d2.RData')


### Analyze KSEI calculations
library(RSiena)
require(network)
require(sna)
require(netseg)
library(Matrix)
load(file='ksei_M230918.RData')
load(file='ksei_M230918_d2.RData')
mydatRF <-  readRDS("rfnetM123_221014.rds") # friendships and race only

# calculate observed and baseline assortativity
# period 1, time 2
fr <- mydatRF$Data1$depvars$friendship[,,2]
ra <- mydatRF$Data1$depvars$raceNet[,,2]
raceCat <- apply(ra, 1, function(z) which(z==1|z==11))
net <- network::network(fr)
network::set.vertex.attribute(net,'race', raceCat)
mm <- mixingmatrix(net,'race')
#netI <- network(1-fr)
#set.vertex.attribute(netI,'race', raceCat)
#mmI <- mixingmatrix(netI,'race') 
ksei_obs <- ei(mm)	
#mm2lay <- array(c(mixingmatrix(netI, 'race'), mm), dim=c(NROW(mm), NCOL(mm), 2))
#or_obs <- orwg(mm2lay)
ksei_base <- rep(0, 1000)
for (i in 1:length(ksei_base)) {
	set.vertex.attribute(net,'race', sample(raceCat, replace=F))
	mm <- mixingmatrix(net,'race')
	ksei_base[i] <- ei(mm)	
	}
# period 2, time 2
fr <- mydatRF$Data2$depvars$friendship[,,2]
ra <- mydatRF$Data2$depvars$raceNet[,,2]
raceCat <- apply(ra, 1, function(z) which(z==1|z==11))
net <- network(fr)
set.vertex.attribute(net,'race', raceCat)
mm <- mixingmatrix(net,'race')
ksei_obs_d2 <- ei(mm)	
ksei_base_d2 <- rep(0, 1000)
for (i in 1:length(ksei_base_d2)) {
	set.vertex.attribute(net,'race', sample(raceCat, replace=F))
	mm <- mixingmatrix(net,'race')
	ksei_base_d2[i] <- ei(mm)	
	}

# create a list for all KSEI objects
ksList <- lapply(ls(pattern="^ksei"), get)
names(ksList) <- ls(pattern="^ksei")
# reverse code the KSEI measure such that -1 = heterophily and +1 = homophily
ksListR <- lapply(ksList, function(x) x*-1)
names(ksListR)
# check that these are equivalent
mean(ksei_inflCre_selAll_d2)
mean(ksList[[14]])
mean(ksList[['ksei_inflCre_selAll_d2']])


# calculate proportions
storeProportions <- matrix(0, nrow=5, ncol=5)
storeProportions[,1] <- c('trend','control','selection','influence','indeterminate')
colnames(storeProportions) <- c('effect','M1','M2','S1','S2')
rangeKSEI <- mean(ksListR[['ksei_full']]) - mean(ksListR[['ksei_base']])
storeProportions[1,2] <- (mean(ksListR[['ksei_trend']]) - mean(ksListR[['ksei_base']])) / rangeKSEI
storeProportions[2,2] <- (mean(ksListR[['ksei_control']]) - mean(ksListR[['ksei_trend']])) / rangeKSEI
storeProportions[3,2] <- min(	(mean(ksListR[['ksei_full']]) - mean(ksListR[['ksei_inflAll_selNon']])),
				(mean(ksListR[['ksei_inflNon_selAll']]) - mean(ksListR[['ksei_control']])) ) / rangeKSEI
storeProportions[4,2] <- min(	(mean(ksListR[['ksei_full']]) - mean(ksListR[['ksei_inflNon_selAll']])),
				(mean(ksListR[['ksei_inflAll_selNon']]) - mean(ksListR[['ksei_control']])) ) / rangeKSEI
storeProportions[5,2] <- 1 - sum(as.numeric(storeProportions[1:4,2]))

rangeKSEI2 <- mean(ksListR[['ksei_full_d2']]) - mean(ksListR[['ksei_base_d2']])
storeProportions[1,3] <- (mean(ksListR[['ksei_trend_d2']]) - mean(ksListR[['ksei_base_d2']])) / rangeKSEI2
storeProportions[2,3] <- (mean(ksListR[['ksei_control_d2']]) - mean(ksListR[['ksei_trend_d2']])) / rangeKSEI2
storeProportions[3,3] <- min(	(mean(ksListR[['ksei_full_d2']]) - mean(ksListR[['ksei_inflAll_selNon_d2']])),
				(mean(ksListR[['ksei_inflNon_selAll_d2']]) - mean(ksListR[['ksei_control_d2']])) ) / rangeKSEI2
storeProportions[4,3] <- min(	(mean(ksListR[['ksei_full_d2']]) - mean(ksListR[['ksei_inflNon_selAll_d2']])),
				(mean(ksListR[['ksei_inflAll_selNon_d2']]) - mean(ksListR[['ksei_control_d2']])) ) / rangeKSEI2
storeProportions[5,3] <- 1 - sum(as.numeric(storeProportions[1:4,3]))

write.table(storeProportions, 'decomp231112.csv', sep=',', row.names=F, quote=F)






##### calculate an OR quickly in igraph/netseg #####
library(igraph)
library(netseg)
#netI <- graph.adjacency(fr, mode='directed')
#netI <- netI %>% set_vertex_attr("race", value = raceCat)
#orwg(netI, 'race')

load('/Users/drschaef/Desktop/workspaceM230927.RData')

# sim data structure
# $sims[[a]][[b]][[c]][[d]]
# a = list, each element is one simulation run (2002) 
# b = list, each element is one period (e.g., Data1, Data2)
# c = list, each element is a dependent network (e.g., friendship, raceNet, ecaEndog)
# d = list, one element, named "1" and contains an edgelist
# begin the decomp
require(network)
require(sna)
require(netseg)
library(Matrix)

# calculate observed and baseline assortativity
# period 1, time 2
fr <- mydatRF$Data1$depvars$friendship[,,2]
ra <- mydatRF$Data1$depvars$raceNet[,,2]
raceCat <- apply(ra, 1, function(z) which(z==1|z==11))
netI <- graph.adjacency(fr, mode='directed')
netI <- netI %>% set_vertex_attr("race", value = raceCat)
or_obs <- orwg(netI, 'race')
pWithinGroup <- sum(table(raceCat) * (table(raceCat)-1)) / (length(raceCat) * (length(raceCat)-1))
(or_base <- pWithinGroup / (1-pWithinGroup))

# a toy model for testing
m_trend <- readRDS('/Users/drschaef/Google Drive/Projects/Adriana/Papers/Homophily/Output/tempM/model230927_trend_M.RDS')
myeff_trendUpdated <- updateTheta(myeff_trend, m_trend)
simAlgorithm30 <- sienaAlgorithmCreate( projname = 'temp', cond=F, useStdInits=F, nsub=0, simOnly=T, n3=30)
s_trend <- siena07( simAlgorithm30, data = mydatRF, effects = myeff_trendUpdated, returnDeps=T)

# recode race into single measure; multiple race get dropped
# but first, count how many have 0 race, more than 1 race
#lapply(simrace20, function(x) table(degree(x, cmode='outdegree')))
extractFrNet <- function(s, DAT) {
	nS <- length(DAT$Data1$nodeSets$students)
	net <- network::network.initialize(n=nS, directed=TRUE)
	edges <- s$'Data1'$'friendship'$'1'[,1:2]
	net <- network::network.edgelist(edges,net)
  	return(net)
	}
extractRaNet <- function(s, DAT) {
	nS <- length(DAT$Data1$nodeSets$students)
	nR <- length(DAT$Data1$nodeSets$erGroups)
	edges <- s$'Data1'$'raceNet'$'1'  
	netDims <- c(nS, nR)  # two mode network dimensions
	edgesPlus <- rbind(edges,c(netDims,0))    # add a final row that ensures the correct network dimensions
	A <- sparseMatrix(i=edgesPlus[,1], j= edgesPlus[,2], x= edgesPlus[,3]) 
	mat <- as.matrix(A)
	return(mat)
	}
calcOR <- function(SIMNET, DAT) {
	simnets <- lapply(SIMNET$sims, function(x) extractFrNet(x,mydatRF))
	simrace <- lapply(SIMNET$sims, function(x) extractRaNet(x,mydatRF))
	OR <- mapply(x=simnets, y=simrace, function(x,y)  {
		nS <- length(DAT$Data1$nodeSets$students)
		od <- sna::degree(y, cmode='outdegree')[1:nS]   # only keep cases reporting one race
		keep <- od==1
		netKeep <- x[keep==1,keep==1]
		raceKeep <- y[keep==1,]
		raceCat <- apply(raceKeep, 1, function(z) which(z==1))
		netI <- graph.adjacency(netKeep, mode='directed')
		netI <- netI %>% set_vertex_attr("race", value = raceCat)
		OR <- orwg(netI, 'race')
		return(OR)		
		})
	return(OR) 
	}

#or_trend <- calcOR(s_trend, mydatRF)  # test
or_trend <- calcOR(m_trend, mydatRF)  # fitted model 
or_control <- calcOR(m_control, mydatRF)  # fitted model
#or_full <- calcOR(s_full, mydatRF)  # first object is full answer object
or_inflAll_selCre <- calcOR(s_inflAll_selCre, mydatRF)  # first object is full answer object
or_inflAll_selEnd <- calcOR(s_inflAll_selEnd, mydatRF)  # first object is full answer object
or_inflAll_selNon <- calcOR(s_inflAll_selNon, mydatRF)  # first object is full answer object
or_inflCre_selAll <- calcOR(s_inflCre_selAll, mydatRF)  # first object is full answer object
or_inflEnd_selAll <- calcOR(s_inflEnd_selAll, mydatRF)  # first object is full answer object
or_inflNon_selAll <- calcOR(s_inflNon_selAll, mydatRF)  # first object is full answer object

setwd('/Users/drschaef/Google Drive/Projects/Adriana/Papers/Homophily/Output')
setwd('/Users/drschaef/My Drive/Projects/Adriana/Papers/Homophily/Output')
#save(list=ls()[ls() %in% ls(pattern="^or")], file='or_M230918.RData')
